options(stringsAsFactors=F)
library(optparse)
library(patchwork)
library(ggplot2)
library(data.table)
library(tidyverse)
library(ggrepel)
library(grid)
library(maftools)
library(ggrepel)
#BiocManager::install("BSgenome.Hsapiens.UCSC.hg19")
library(BSgenome.Hsapiens.UCSC.hg19)

##########################################################################################

option_list <- list(
  make_option(c("--gisticAllLesionsFile_IGC"), type = "character"),
  make_option(c("--gisticAmpGenesFile_IGC"), type = "character"),
  make_option(c("--gisticDelGenesFile_IGC"), type = "character"),
  make_option(c("--gisticScoresFile_IGC"), type = "character"),
  make_option(c("--gisticAllLesionsFile_DGC"), type = "character"),
  make_option(c("--gisticAmpGenesFile_DGC"), type = "character"),
  make_option(c("--gisticDelGenesFile_DGC"), type = "character"),
  make_option(c("--gisticScoresFile_DGC"), type = "character"),
  make_option(c("--show_Cytoband_amp_file"), type = "character"),
  make_option(c("--show_Cytoband_del_file"), type = "character"),
  make_option(c("--out_path"), type = "character"), 
  make_option(c("--Class"), type = "character")
)

if(1!=1){
  
  gisticAllLesionsFile_IGC = 'all_lesions.conf_99.txt'
  gisticAmpGenesFile_IGC = 'amp_genes.conf_99.txt'
  gisticDelGenesFile_IGC = 'del_genes.conf_99.txt'
  gisticScoresFile_IGC = 'scores.gistic_IGC'
  
}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gisticAllLesionsFile_IGC <- opt$gisticAllLesionsFile_IGC
gisticAmpGenesFile_IGC <- opt$gisticAmpGenesFile_IGC
gisticDelGenesFile_IGC <- opt$gisticDelGenesFile_IGC
gisticScoresFile_IGC <- opt$gisticScoresFile_IGC
gisticAllLesionsFile_DGC <- opt$gisticAllLesionsFile_DGC
gisticAmpGenesFile_DGC <- opt$gisticAmpGenesFile_DGC
gisticDelGenesFile_DGC <- opt$gisticDelGenesFile_DGC
gisticScoresFile_DGC <- opt$gisticScoresFile_DGC
show_Cytoband_amp_file <- opt$show_Cytoband_amp_file
show_Cytoband_del_file <- opt$show_Cytoband_del_file
out_path <- opt$out_path
Class <- opt$Class

###########################################################################################
#建立df，包含各染色体的长度信息

df <- data.frame(chromName = seqnames(BSgenome.Hsapiens.UCSC.hg19), # 染色体名字
                 chromlength = seqlengths(BSgenome.Hsapiens.UCSC.hg19)# 每条染色体长度
)

df$chromNum <- 1:length(df$chromName) # 染色体名字纯数字版，为了和scores.gistic里面的对应
df <- df[1:22,] # 只要前22条染色体信息


df$chromlengthCumsum <- cumsum(as.numeric(df$chromlength)) # 染色体累加长度

# 得到每条染色体从0开始的起始坐标
df$chormStartPosFrom0 <- c(0,df$chromlengthCumsum[-nrow(df)])

# 计算每条染色体中间位置坐标，用来最后加文字
tmp_middle <- diff(c(0,df$chromlengthCumsum)) / 2
df$chromMidelePosFrom0 <- df$chormStartPosFrom0 + tmp_middle

show_Cytoband_amp <- unique(data.frame(fread(show_Cytoband_amp_file))$Cytoband)
show_Cytoband_del <- unique(data.frame(fread(show_Cytoband_del_file))$Cytoband)
show_Cytoband <- unique(c(show_Cytoband_amp , show_Cytoband_del))

###########################################################################################
##设置的函数

gisticChromGet = function(gistic = gistic, fdrCutOff = fdrCutOff, ref.build = ref.build,class=class,sub=sub) {
  
g = maftools::getCytobandSummary(gistic)
g = g[qvalues < fdrCutOff]
g[,Chromosome := sapply(strsplit(x = g$Wide_Peak_Limits, split = ':'), '[', 1)]
g[,loc := sapply(strsplit(x = g$Wide_Peak_Limits, split = ':'), '[', 2)]
g[,Start_Position := sapply(strsplit(x = g$loc, split = '-'), '[', 1)]
g[,End_Position := sapply(strsplit(x = g$loc, split = '-'), '[', 2)]
g.lin = maftools:::transformSegments(segmentedData = g[,.(Chromosome, Start_Position, End_Position, qvalues, Cytoband, Variant_Classification)])

g.lin <- subset(g.lin,Variant_Classification==sub)
mb <- g.lin
#mb$StartPos <- mb$Start_Position_updated
#mb$EndPos <- mb$End_Position_updated

gis.scores = maftools:::transformSegments(segmentedData = gistic@gis.scores, build = ref.build)
gis.scores <- subset(gis.scores,Variant_Classification==sub)
#gis.scores$amp = ifelse(test = gis.scores$Variant_Classification == 'Del', yes = -gis.scores$fdr, no = gis.scores$fdr)
gis.scores$amp = gis.scores$G_Score

fdrCutOff = -log10(fdrCutOff)

##设置主键
data.table::setkey(x = gis.scores, Chromosome, Start_Position_updated, End_Position_updated)
cyto_peaks_scores = data.table::foverlaps(y = gis.scores[,.(Chromosome, Start_Position_updated, End_Position_updated, amp,G_Score,fdr)],
                                          x = mb[,.(Chromosome, Start_Position_updated, End_Position_updated, Cytoband,qvalues)])
cyto_peaks_scores = cyto_peaks_scores[order(Cytoband, abs(amp), decreasing = TRUE)][Cytoband %in% mb$Cytoband][!duplicated(Cytoband)]
cyto_peaks_scores = cyto_peaks_scores[complete.cases(cyto_peaks_scores)]

gis.scores$name <- paste0(gis.scores$Chromosome,"_",gis.scores$Start_Position_updated,"_",gis.scores$End_Position_updated)
cyto_peaks_scores$name <- paste0(cyto_peaks_scores$Chromosome,"_",cyto_peaks_scores$Start_Position_updated,"_",cyto_peaks_scores$End_Position_updated)
gis_all <- merge(gis.scores,cyto_peaks_scores,by="name",all.x=T)
gis_all$class <- class
return(gis_all)

}

plotChrom <- function(df=df,max_x=max_x){
  p_rect=ggplot()+
    geom_rect(data = df[df$chromNum %% 2 == 0,],
              aes(
                xmin = chormStartPosFrom0,
                xmax = chromlengthCumsum,
                ymin = -0.01,ymax = 0
              ),
              inherit.aes = FALSE,
              col="white",
              fill = "white",
              show.legend = F)+  
    geom_rect(data = df[df$chromNum %% 2 != 0,],
              aes(
                xmin = chormStartPosFrom0,
                xmax = chromlengthCumsum,
                ymin = -0.01,ymax = 0
              ),
              inherit.aes = FALSE,
              col="#D3D3D3",
              fill="#D3D3D3",
              show.legend = F)+
    geom_text(data = df,aes(x=chromMidelePosFrom0,y=-0.005,label=chromNum),col="black")+
    scale_y_continuous(limits = c(-0.01,0),labels = NULL)+
    scale_x_continuous(limits = c(0, max_x), expand = c(0, 0)) +
    theme_minimal()+
    ylab("")+
    theme( panel.grid.major = element_blank(),
           panel.grid.minor = element_blank(),
           panel.background = element_blank(),
           axis.text.x = element_blank(),
           axis.title.x = element_blank(),
           axis.ticks.x = element_blank(),
           axis.text.y = element_text(color = "black", size = 12),
           axis.title = element_text(color = "black", size = 12),
           legend.title = element_blank(),
           legend.text = element_text(color = "black", size = 10)) 
}


scale_y_continuous_set <- function(uplimit=uplimit,class=class,sep=sep) {
  if (class == "IGC") {
    return(scale_y_continuous(
      expand = c(0.02, 0.02),
      limits = c(0, uplimit),
      breaks = seq(0, uplimit, by = sep)
    ))
  } else if(class == "DGC"){
    return(scale_y_continuous(
      expand = c(0.02, 0.02),
      limits = c(-uplimit, 0),
      breaks = seq(0,-uplimit, by = -sep)
    ))
  }
}

plotFocal <- function(gistic=gistic,up_limit=up_limit,down_limit=down_limit,Cytoband_limit=Cytoband_limit,color=color,max_x=max_x,class=class,sep=sep){
  gistic$Cytoband <- ifelse(gistic$Cytoband %in% show_Cytoband , gistic$Cytoband , "")
  p <- ggplot(gistic, aes(Start_Position_updated.x, amp.x)) +
    geom_rect(data = df[df$chromNum %% 2 == 0,],
              aes(
                xmin = chormStartPosFrom0,
                xmax = chromlengthCumsum,
                ymin = -Inf,ymax = Inf
              ),
              inherit.aes = FALSE,
              col="white",
              fill = "white",
              show.legend = F)+  
    geom_rect(data = df[df$chromNum %% 2 != 0,],
              aes(
                xmin = chormStartPosFrom0,
                xmax = chromlengthCumsum,
                ymin = -Inf,ymax = Inf
              ),
              inherit.aes = FALSE,
              col="#D3D3D3",
              fill="#D3D3D3",
              show.legend = F)+
    geom_area(aes(group=class,fill=class,colour = class)) +
    scale_colour_manual(values = color ) +
    scale_fill_manual(values=scales::alpha("white", 0))+
    geom_text_repel(data = subset(gistic, !is.na(Cytoband) & abs(amp.x) > Cytoband_limit & abs(qvalues) < 0.2),aes(x = i.Start_Position_updated, y = amp.x, label = Cytoband),force = 0.1, color = "grey20",size = 2.5, point.padding = 0.5, arrow = arrow(length = unit(0.02, "npc"), type = "open", ends = "last"), max.iter = 1000)+
    geom_vline(data = df ,mapping=aes(xintercept=chromlengthCumsum),linetype=2)+
    #geom_hline(yintercept = log10(0.25), linetype = 2)+
    #geom_hline(yintercept = -log10(0.25), linetype = 2)+
    ylab("G-score")+
	scale_y_continuous_set(uplimit=up_limit,class=class,sep=sep)+
    scale_x_continuous(limits = c(0, max_x), expand = c(0, 0)) +
    theme_classic()+
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_blank(),
          axis.text.x = element_blank(),
          axis.title.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.text.y = element_text(color = "black", size = 12),
          axis.title = element_text(color = "black", size = 12),
          legend.title = element_blank(),
          legend.text = element_text(color = "black", size = 10))
}

###########################################################################################
##导入数据
gistic_IGC = readGistic(gisticAllLesionsFile = gisticAllLesionsFile_IGC,
                        gisticAmpGenesFile = gisticAmpGenesFile_IGC, 
                        gisticDelGenesFile = gisticDelGenesFile_IGC, 
                        gisticScoresFile = gisticScoresFile_IGC, 
                        isTCGA =F)

gistic_DGC = readGistic(gisticAllLesionsFile = gisticAllLesionsFile_DGC,
                        gisticAmpGenesFile = gisticAmpGenesFile_DGC, 
                        gisticDelGenesFile = gisticDelGenesFile_DGC, 
                        gisticScoresFile = gisticScoresFile_DGC, 
                        isTCGA =F)

###########################################################################################
##数据整合(amp)

gis_all_IGC <- gisticChromGet(gistic = gistic_IGC, fdrCutOff = 0.25, ref.build = "hg19",class="IGC",sub="Amp")
gis_all_DGC <- gisticChromGet(gistic = gistic_DGC, fdrCutOff = 0.25, ref.build = "hg19",class="DGC",sub="Amp")

gis_all_DGC$amp.x <- -gis_all_DGC$amp.x
gis_all_DGC$amp.y <- -gis_all_DGC$amp.y


###########################################################################################
##画图（amp）
color <- c("#A32A31","#1e69b4")
names(color) <- c("IGC" , "DGC")

up_limit <- max(gis_all_IGC$amp.x)+0.5
down_limit <- min(gis_all_DGC$amp.x)-0.5
max_x <- max(df$chromlengthCumsum)
print(down_limit)

p1 <- plotFocal(gistic=gis_all_IGC,up_limit=up_limit,max_x=max_x,Cytoband_limit=0.5, color=color,class="IGC",sep=0.5)
p2 <- plotFocal(gistic=gis_all_DGC,up_limit=up_limit,max_x=max_x,Cytoband_limit=0.5, color=color,class="DGC",sep=0.5)
p_rect <- plotChrom(df=df,max_x=max_x)
p <- p1/p_rect/p2+plot_layout(heights = c(1, 0.2, 1))

out_name <- paste0(out_path,"/amp/CIN_IGCcompareDGC_focal_amp_",Class,".gscore.pdf")
ggsave(out_name, plot = p, width = 15, height = 6)


###########################################################################################
##数据整合(del)

gis_all_IGC_del <- gisticChromGet(gistic = gistic_IGC, fdrCutOff = 0.25, ref.build = "hg19",class="IGC",sub="Del")
gis_all_DGC_del <- gisticChromGet(gistic = gistic_DGC, fdrCutOff = 0.25, ref.build = "hg19",class="DGC",sub="Del")

gis_all_DGC_del$amp.x <- -gis_all_DGC_del$amp.x
gis_all_DGC_del$amp.y <- -gis_all_DGC_del$amp.y

write.table(gis_all_IGC_del,paste0(out_path,"/del/test.tsv"),row.names=F,sep="\t",quote=F)
write.table(gis_all_DGC_del,paste0(out_path,"/del/test_2.tsv"),row.names=F,sep="\t",quote=F)
###########################################################################################
##画图（Del）
color <- c("#A32A31","#1e69b4")
names(color) <- c("IGC" , "DGC")

up_limit <- max(gis_all_IGC_del$amp.x)+0.3
down_limit <- min(gis_all_DGC_del$amp.x)-0.3
max_x <- max(df$chromlengthCumsum)

#print(subset(gis_all_DGC_del, !is.na(Cytoband) & abs(amp.x) > 4))

p1 <- plotFocal(gistic=gis_all_IGC_del,up_limit=up_limit,max_x=max_x,Cytoband_limit=0.5, color=color,class="IGC",sep=0.3)
p2 <- plotFocal(gistic=gis_all_DGC_del,up_limit=up_limit,max_x=max_x,Cytoband_limit=0.5, color=color,class="DGC",sep=0.3)
p_rect <- plotChrom(df=df,max_x=max_x)
p <- p1/p_rect/p2+plot_layout(heights = c(1, 0.2, 1))

out_name <- paste0(out_path,"/del/CIN_IGCcompareDGC_focal_del_",Class,".gscore.pdf")
ggsave(out_name, plot = p, width = 15, height = 6)

